Deep dive regional communities - using relative abundance prox

source('../env.R')

Species in communities

It seems reasonable to expect that cities with simialr regional pools will have similar species entering the city, and thus a similar response to urbanisation.

Load data

communities = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'communities_for_analysis.csv'))

Rows: 2462 Columns: 10── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): city_name, jetz_species_name, seasonal, presence, origin
dbl (2): city_id, relative_abundance_proxy
lgl (3): present_urban_high, present_urban_med, present_urban_low
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

communities_summary = communities %>% group_by(city_id) %>% summarise(
  regional_pool_size = n(), 
  urban_pool_size = sum(relative_abundance_proxy > 0)
)

ggplot(communities %>% filter(relative_abundance_proxy > 0), aes(x = relative_abundance_proxy)) + geom_bar(stat = "bin")

city_points = st_centroid(read_sf(filename(CITY_DATA_OUTPUT_DIR, 'city_selection.shp')))

Warning: st_centroid assumes attributes are constant over geometriesWarning: st_centroid does not give correct centroids for longitude/latitude data

community_data_metrics = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'community_assembly_metrics_using_relative_abundance.csv')) %>%
  dplyr::select(city_id, mntd_normalised, fdiv_normalised, mass_fdiv_normalised, locomotory_trait_fdiv_normalised, trophic_trait_fdiv_normalised, gape_width_fdiv_normalised) %>%
  left_join(read_csv(filename(CITY_DATA_OUTPUT_DIR, 'realms.csv'))) %>%
  left_join(communities_summary) %>%
  left_join(city_points[,c('city_id', 'city_nm')]) %>%
  rename(city_name='city_nm') %>%
  na.omit() %>%
  arrange(city_id)

Rows: 341 Columns: 37── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
dbl (37): mntd_normalised, mntd_actual, mntd_min, mntd_max, mntd_mean, mntd_sd, fdiv_normalised, fdiv_actual, fdiv_min, fdiv_max, fdiv_...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Rows: 342 Columns: 2── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): core_realm
dbl (1): city_id
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Joining with `by = join_by(city_id)`Joining with `by = join_by(city_id)`Joining with `by = join_by(city_id)`

community_data_metrics

Load trait data

traits = read_csv(filename(TAXONOMY_OUTPUT_DIR, 'traits_jetz.csv'))

Rows: 304 Columns: 5── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): jetz_species_name
dbl (4): gape_width, trophic_trait, locomotory_trait, mass
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(traits)

fetch_normalised_traits = function(required_species_list) {
  required_traits = traits %>% filter(jetz_species_name %in% required_species_list)
  
  required_traits$gape_width_normalised = normalise(required_traits$gape_width, min(required_traits$gape_width), max(required_traits$gape_width))
  required_traits$trophic_trait_normalised = normalise(required_traits$trophic_trait, min(required_traits$trophic_trait), max(required_traits$trophic_trait))
  required_traits$locomotory_trait_normalised = normalise(required_traits$locomotory_trait, min(required_traits$locomotory_trait), max(required_traits$locomotory_trait))
  required_traits$mass_normalised = normalise(required_traits$mass, min(required_traits$mass), max(required_traits$mass))
  
  traits_normalised_long = required_traits %>% pivot_longer(cols = c('gape_width_normalised', 'trophic_trait_normalised', 'locomotory_trait_normalised', 'mass_normalised'), names_to = 'trait', values_to = 'normalised_value') %>% dplyr::select(jetz_species_name, trait, normalised_value)
  traits_normalised_long$trait = factor(traits_normalised_long$trait, levels = c('gape_width_normalised', 'trophic_trait_normalised', 'locomotory_trait_normalised', 'mass_normalised'), labels = c('Gape Width', 'Trophic Trait', 'Locomotory Trait', 'Mass'))
  
  traits_normalised_long
}

fetch_normalised_traits(c('Aplopelia_larvata', 'Chalcophaps_indica', 'Caloenas_nicobarica'))

Read in our phylogeny

phylo_tree = read.tree(filename(TAXONOMY_OUTPUT_DIR, 'phylogeny.tre'))
ggtree(phylo_tree, layout='circular')

Load resolve ecoregions

resolve = read_resolve()

Warning: st_buffer does not correctly buffer longitude/latitude datadist is assumed to be in decimal degrees (arc_degrees).
Warning: st_simplify does not correctly simplify longitude/latitude data, dTolerance needs to be in decimal degrees

Create helper functions

to_species_matrix = function(filtered_communities) {
  filtered_communities %>% 
    dplyr::select(city_id, jetz_species_name) %>% 
    distinct() %>%
    mutate(present = TRUE) %>% 
    pivot_wider(
      names_from = jetz_species_name, 
      values_from = "present", 
      values_fill = list(present = F)
    ) %>% 
    tibble::column_to_rownames(var='city_id')
}

community_nmds = function(filtered_communities) {
  species_matrix = to_species_matrix(filtered_communities)
  nmds <- metaMDS(species_matrix, k=2, trymax = 30)
  nmds_result = data.frame(scores(nmds)$sites)
  nmds_result$city_id = as.double(rownames(nmds_result))
  rownames(nmds_result) = NULL
  nmds_result
}

https://www.datacamp.com/tutorial/k-means-clustering-r

scree_plot = function(community_nmds_data) {
  # Decide how many clusters to look at
  n_clusters <- min(10, nrow(community_nmds_data) - 1)
  
  # Initialize total within sum of squares error: wss
  wss <- numeric(n_clusters)
  
  set.seed(123)
  
  # Look over 1 to n possible clusters
  for (i in 1:n_clusters) {
    # Fit the model: km.out
    km.out <- kmeans(community_nmds_data[,c('NMDS1','NMDS2')], centers = i, nstart = 20)
    # Save the within cluster sum of squares
    wss[i] <- km.out$tot.withinss
  }
  
  # Produce a scree plot
  wss_df <- tibble(clusters = 1:n_clusters, wss = wss)
   
  scree_plot <- ggplot(wss_df, aes(x = clusters, y = wss, group = 1)) +
      geom_point(size = 4) +
      geom_line() +
      geom_hline(linetype="dashed", color = "orange", yintercept = wss) +
      scale_x_continuous(breaks = c(2, 4, 6, 8, 10)) +
      xlab('Number of clusters')
  scree_plot
}

cluster_cities = function(city_nmds, cities_community_data, centers) {
  set.seed(123)
  kmeans_clusters <- kmeans(city_nmds[,c('NMDS1', 'NMDS2')], centers = centers, nstart = 20)
  city_nmds$cluster = kmeans_clusters$cluster
  cities_community_data %>% left_join(city_nmds) %>% mutate(cluster = as.factor(cluster))
}

plot_nmds_clusters = function(cluster_cities) {
  cluster_cities %>% dplyr::select(city_id, city_name, NMDS1, NMDS2, cluster) %>% distinct() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, colour = cluster)) + geom_point() + geom_label_repel(aes(label = city_name))
}

plot_city_cluster = function(city_cluster_data_metrics, title) {
  species_in_cluster = communities %>% 
    filter(city_id %in% city_cluster_data_metrics$city_id) %>% 
    dplyr::select(jetz_species_name, city_name, relative_abundance_proxy)

  
  tree_cropped <- ladderize(drop.tip(phylo_tree, setdiff(phylo_tree$tip.label, species_in_cluster$jetz_species_name)))
    
  gg_tree = ggtree(tree_cropped)
  
  gg_presence = ggplot(species_in_cluster, aes(x=city_name, y=jetz_species_name)) + 
          geom_tile(aes(fill=relative_abundance_proxy)) + 
          scale_fill_gradientn(colours=c("#98FB98", "#FFFFE0", "yellow", "orange", "#FF4500", "red", "red"), values=c(0, 0.00000000001, 0.1, 0.25, 0.5, 0.75, 1), na.value = "transparent") +
          theme_minimal() + xlab(NULL) + ylab(NULL) + 
          theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
          labs(fill='Urban Proxy Abundance')
  
  species_in_cluster_traits = fetch_normalised_traits(species_in_cluster$jetz_species_name)
  
  gg_traits = ggplot(species_in_cluster_traits, aes(x = trait, y = jetz_species_name, size = normalised_value)) + geom_point() + theme_minimal() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), axis.text.y=element_blank()) + xlab(NULL) + ylab(NULL) + labs(size = "Normalised Value")
  
  gg_cities_mntd = ggplot(city_cluster_data_metrics, aes(x = city_name, y = mntd_normalised)) + geom_bar(stat = "identity") + theme_minimal() + theme(legend.position = "none", axis.text.x=element_blank()) + xlab(NULL) + ylab("MNTD") + ylim(0, 1)
  
  gg_cities_fd = ggplot(city_cluster_data_metrics, aes(x = city_name, y = fdiv_normalised)) + geom_bar(stat = "identity") + theme_minimal() + theme(legend.position = "none", axis.text.x=element_blank()) + xlab(NULL) + ylab("FDiv") + ylim(0, 1)
  
  gg_title = ggplot() + labs(title = title) + theme_minimal()
  
  gg_presence %>% insert_top(gg_cities_mntd, height = 0.5) %>% insert_top(gg_cities_fd, height = 0.5) %>% insert_left(gg_tree, width = 0.75) %>% insert_right(gg_traits, width = 0.5) %>% insert_top(gg_title, height = 0.06)
}

REGION_DEEP_DIVE_FIGURES_OUTPUT = mkdir(FIGURES_OUTPUT_DIR, 'appendix_regional_deep_dive_using_abundance')

Nearctic

nearctic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Nearctic')
nearctic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()

$city_name
 [1] "San Jose"                 "Los Angeles"              "Concord"                  "Tijuana"                  "Bakersfield"             
 [6] "Fresno"                   "Sacramento"               "Mexicali"                 "Hermosillo"               "Las Vegas"               
[11] "Phoenix"                  "Tucson"                   "Durango"                  "Portland"                 "Chihuahua"               
[16] "Aguascalientes"           "Seattle"                  "Ciudad Juárez"            "San Luis Potosí"          "Mexico City"             
[21] "Saltillo"                 "Vancouver"                "Salt Lake City"           "Albuquerque"              "Monterrey"               
[26] "Nuevo Laredo"             "San Antonio"              "Denver"                   "Austin"                   "Houston"                 
[31] "Dallas"                   "Oklahoma City"            "Calgary"                  "New Orleans"              "Kansas City"             
[36] "Omaha"                    "St. Louis"                "Bradenton"                "Tampa"                    "Minneapolis [Saint Paul]"
[41] "Atlanta"                  "Orlando"                  "Louisville"               "Chicago"                  "Indianapolis"            
[46] "Milwaukee"               

attr(,"na.action")
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
attr(,"class")
[1] "omit"

nearctic_cities_nmds = community_nmds(communities %>% filter(city_id %in% nearctic_cities_community_data$city_id)) 

Run 0 stress 0.1005678 
Run 1 stress 0.1210166 
Run 2 stress 0.1012776 
Run 3 stress 0.1007176 
... Procrustes: rmse 0.009402672  max resid 0.03496107 
Run 4 stress 0.1000046 
... New best solution
... Procrustes: rmse 0.007232701  max resid 0.0346979 
Run 5 stress 0.1228731 
Run 6 stress 0.1022505 
Run 7 stress 0.1205421 
Run 8 stress 0.1205421 
Run 9 stress 0.1240951 
Run 10 stress 0.1217438 
Run 11 stress 0.1000046 
... New best solution
... Procrustes: rmse 0.000004447461  max resid 0.00001302201 
... Similar to previous best
Run 12 stress 0.1000046 
... Procrustes: rmse 0.000003090941  max resid 0.00001098798 
... Similar to previous best
Run 13 stress 0.121519 
Run 14 stress 0.1012776 
Run 15 stress 0.1234703 
Run 16 stress 0.1217145 
Run 17 stress 0.1007176 
Run 18 stress 0.1217145 
Run 19 stress 0.1205423 
Run 20 stress 0.124095 
*** Best solution repeated 2 times

nearctic_cities_nmds

scree_plot(nearctic_cities_nmds)

nearctic_cities = cluster_cities(city_nmds = nearctic_cities_nmds, cities_community_data = nearctic_cities_community_data, centers = 5)

Joining with `by = join_by(city_id)`

plot_nmds_clusters(nearctic_cities)

nearctic_biomes = st_crop(resolve[resolve$REALM == 'Nearctic',c('REALM')], xmin = -220, ymin = 0, xmax = 0, ymax = 70)

although coordinates are longitude/latitude, st_intersection assumes that they are planar
Warning: attribute variables are assumed to be spatially constant throughout all geometries

 
ggplot() + 
  geom_sf(data = nearctic_biomes, aes(geometry = geometry)) + 
  geom_sf(data = nearctic_cities, aes(geometry = geometry, color = cluster))

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_clusters.jpg'))

Saving 7.29 x 4.51 in image

Neartic Cluster 1`

nearctic_cities %>% filter(cluster == 1) %>% plot_city_cluster('Neartic cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster1.jpg'))

Saving 7.29 x 4.51 in image

Neartic Cluster 2

nearctic_cities %>% filter(cluster == 2) %>% plot_city_cluster('Neartic cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster2.jpg'))

Saving 7.29 x 4.51 in image

Neartic Cluster 3

nearctic_cities %>% filter(cluster == 3) %>% plot_city_cluster('Neartic cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster3.jpg'))

Saving 7.29 x 4.51 in image

Neartic Cluster 4

nearctic_cities %>% filter(cluster == 4) %>% plot_city_cluster('Neartic cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster4.jpg'))

Saving 7.29 x 4.51 in image

Neartic Cluster 5

nearctic_cities %>% filter(cluster == 5) %>% plot_city_cluster('Neartic cluster 5')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster5.jpg'))

Saving 7.29 x 4.51 in image

Neotropic

neotropic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Neotropic')
neotropic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()

$city_name
 [1] "Culiacán"                  "Guadalajara"               "Morelia"                   "Acapulco"                 
 [5] "Querétaro"                 "Cuernavaca"                "Puebla"                    "Oaxaca"                   
 [9] "Xalapa"                    "Veracruz"                  "Tuxtla Gutiérrez"          "Villahermosa"             
[13] "Guatemala City"            "San Salvador"              "San Pedro Sula"            "Mérida"                   
[17] "Tegucigalpa"               "Managua"                   "San José"                  "Cancún"                   
[21] "Guayaquil"                 "Chiclayo"                  "Panama City"               "Trujillo"                 
[25] "Quito"                     "Havana"                    "Cali"                      "Lima"                     
[29] "Pereira"                   "Miami"                     "Medellín"                  "Ibagué"                   
[33] "Cartagena"                 "Kingston"                  "Bogota"                    "Barranquilla"             
[37] "Bucaramanga"               "Cúcuta"                    "Maracaibo"                 "Arequipa"                 
[41] "Barquisimeto"              "Santo Domingo"             "Maracay"                   "El Alto [La Paz]"         
[45] "Caracas"                   "Cochabamba"                "Viña del Mar [Valparaíso]" "Río Piedras [San Juan]"   
[49] "Barcelona"                 "Concepción"                "Santiago"                  "Mendoza"                  
[53] "Salta"                     "Cordoba"                   "Asuncion"                  "Buenos Aires"             
[57] "La Plata"                  "Ciudad del Este"           "Montevideo"                "Mar del Plata"            
[61] "Porto Alegre"              "São Paulo"                 "Santos"                    "Sao Jose dos Campos"      

attr(,"na.action")
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
attr(,"class")
[1] "omit"

neotropic_cities_nmds = community_nmds(communities %>% filter(city_id %in% neotropic_cities_community_data$city_id)) 

Run 0 stress 0.134619 
Run 1 stress 0.1414255 
Run 2 stress 0.1348046 
... Procrustes: rmse 0.0106512  max resid 0.05474793 
Run 3 stress 0.1348237 
... Procrustes: rmse 0.01071367  max resid 0.05473745 
Run 4 stress 0.1405639 
Run 5 stress 0.1414222 
Run 6 stress 0.141222 
Run 7 stress 0.1414222 
Run 8 stress 0.1344509 
... New best solution
... Procrustes: rmse 0.006925358  max resid 0.04568131 
Run 9 stress 0.1402357 
Run 10 stress 0.1348046 
... Procrustes: rmse 0.006662375  max resid 0.04608377 
Run 11 stress 0.1405635 
Run 12 stress 0.1405642 
Run 13 stress 0.134619 
... Procrustes: rmse 0.006938373  max resid 0.04565142 
Run 14 stress 0.1346368 
... Procrustes: rmse 0.006829873  max resid 0.04569156 
Run 15 stress 0.1405644 
Run 16 stress 0.1348237 
... Procrustes: rmse 0.006540716  max resid 0.04604158 
Run 17 stress 0.134433 
... New best solution
... Procrustes: rmse 0.001186138  max resid 0.006984211 
... Similar to previous best
Run 18 stress 0.1405635 
Run 19 stress 0.1344509 
... Procrustes: rmse 0.001186071  max resid 0.006831159 
... Similar to previous best
Run 20 stress 0.1344509 
... Procrustes: rmse 0.001186504  max resid 0.006792118 
... Similar to previous best
*** Best solution repeated 3 times

neotropic_cities_nmds

scree_plot(neotropic_cities_nmds)

neotropic_cities = cluster_cities(city_nmds = neotropic_cities_nmds, cities_community_data = neotropic_cities_community_data, centers = 5)

Joining with `by = join_by(city_id)`

plot_nmds_clusters(neotropic_cities)

neotropic_biomes = resolve[resolve$REALM == 'Neotropic',c('REALM')]
 
ggplot() + 
  geom_sf(data = neotropic_biomes, aes(geometry = geometry)) + 
  geom_sf(data = neotropic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_clusters.jpg'))

Saving 7.29 x 4.51 in image

Neotropic Cluster 1

neotropic_cities %>% filter(cluster == 1) %>% plot_city_cluster('Neotropic cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster1.jpg'))

Saving 7.29 x 4.51 in image

Neotropic Cluster 2

neotropic_cities %>% filter(cluster == 2) %>% plot_city_cluster('Neotropic cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster2.jpg'))

Saving 7.29 x 4.51 in image

Neotropic Cluster 3

neotropic_cities %>% filter(cluster == 3) %>% plot_city_cluster('Neotropic cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster3.jpg'))

Saving 7.29 x 4.51 in image

Neotropic Cluster 4

neotropic_cities %>% filter(cluster == 4) %>% plot_city_cluster('Neotropic cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster4.jpg'))

Saving 7.29 x 4.51 in image

Neotropic Cluster 5

neotropic_cities %>% filter(cluster == 5) %>% plot_city_cluster('Neotropic cluster 5')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster5.jpg'))

Saving 7.29 x 4.51 in image

Palearctic

palearctic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Palearctic')
palearctic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()

$city_name
 [1] "Lisbon"                "Porto"                 "Marrakesh"             "Seville"               "Dublin"               
 [6] "Málaga"                "Madrid"                "Glasgow"               "Bilbao"                "Liverpool"            
[11] "Bristol"               "Manchester"            "Birmingham"            "Leeds"                 "Newcastle upon Tyne"  
[16] "Sheffield"             "Nottingham"            "Valencia"              "London"                "Toulouse"             
[21] "Paris"                 "Barcelona"             "Rotterdam [The Hague]" "Brussels"              "Amsterdam"            
[26] "Lyon"                  "Marseille"             "Dusseldorf"            "Nice"                  "Frankfurt am Main"    
[31] "Zurich"                "Oslo"                  "Stuttgart"             "Hamburg"               "Genoa"                
[36] "Nuremberg"             "Copenhagen"            "Munich"                "Berlin"                "Dresden"              
[41] "Rome"                  "Prague"                "Stockholm"             "Poznan"                "Vienna"               
[46] "Wroclaw"               "Zagreb"                "Gdansk"                "Budapest"              "Krakow"               
[51] "Warsaw"                "Helsinki"              "Riga"                  "Belgrade"              "Lviv"                 
[56] "Sofia"                 "Thessaloniki"          "Saint Petersburg"      "Minsk"                 "Athens"               
[61] "Kyiv"                  "Istanbul"              "Odesa"                 "Samsun"                "Luxor"                
[66] "Tel Aviv"              "Jerusalem"             "Tbilisi"               "Yerevan"               "Kuwait City"          
[71] "Doha"                  "Abu Dhabi"             "Dubai"                 "Bishkek"              

attr(,"na.action")
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
attr(,"class")
[1] "omit"

palearctic_cities_nmds = community_nmds(communities %>% filter(city_id %in% palearctic_cities_community_data$city_id)) 

Run 0 stress 0.04961857 
Run 1 stress 0.06079762 
Run 2 stress 0.06489516 
Run 3 stress 0.05546046 
Run 4 stress 0.06348958 
Run 5 stress 0.05395553 
Run 6 stress 0.04967654 
... Procrustes: rmse 0.02445083  max resid 0.1151427 
Run 7 stress 0.06304717 
Run 8 stress 0.07947601 
Run 9 stress 0.05001033 
... Procrustes: rmse 0.06742835  max resid 0.2224835 
Run 10 stress 0.05392775 
Run 11 stress 0.06304664 
Run 12 stress 0.05258764 
Run 13 stress 0.08298216 
Run 14 stress 0.05660148 
Run 15 stress 0.05105666 
Run 16 stress 0.05672742 
Run 17 stress 0.06304727 
Run 18 stress 0.05479118 
Run 19 stress 0.05101845 
Run 20 stress 0.05471648 
Run 21 stress 0.07604128 
Run 22 stress 0.08017294 
Run 23 stress 0.06304694 
Run 24 stress 0.04961852 
... New best solution
... Procrustes: rmse 0.0000829673  max resid 0.0003983299 
... Similar to previous best
*** Best solution repeated 1 times

palearctic_cities_nmds

scree_plot(palearctic_cities_nmds)

palearctic_cities = cluster_cities(city_nmds = palearctic_cities_nmds, cities_community_data = palearctic_cities_community_data, centers = 6)

Joining with `by = join_by(city_id)`

plot_nmds_clusters(palearctic_cities)

palearctic_biomes = resolve[resolve$REALM == 'Palearctic',c('REALM')]
 
ggplot() + 
  geom_sf(data = palearctic_biomes, aes(geometry = geometry)) + 
  geom_sf(data = palearctic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_clusters.jpg'))

Saving 7.29 x 4.51 in image

Palearctic Cluster 1

palearctic_cities %>% filter(cluster == 1) %>% plot_city_cluster('Palearctic cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster1.jpg'))

Saving 7.29 x 4.51 in image

Palearctic Cluster 2

palearctic_cities %>% filter(cluster == 2) %>% plot_city_cluster('Palearctic cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster2.jpg'))

Saving 7.29 x 4.51 in image

Palearctic Cluster 3

palearctic_cities %>% filter(cluster == 3) %>% plot_city_cluster('Palearctic cluster 3')

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster3.jpg'))

Saving 14 x 6 in image

Palearctic Cluster 4

palearctic_cities %>% filter(cluster == 4) %>% plot_city_cluster('Palearctic cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster4.jpg'))

Saving 7.29 x 4.51 in image

Palearctic Cluster 5

palearctic_cities %>% filter(cluster == 5) %>% plot_city_cluster('Palearctic cluster 5')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster5.jpg'))

Saving 7.29 x 4.51 in image

Palearctic Cluster 6

palearctic_cities %>% filter(cluster == 6) %>% plot_city_cluster('Palearctic cluster 6')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster6.jpg'))

Saving 7.29 x 4.51 in image

Afrotropic

afrotropic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Afrotropic')
afrotropic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()

$city_name
[1] "Cape Town"    "Johannesburg" "Pretoria"     "Kigali"       "Kampala"      "Arusha"       "Nairobi"      "Addis Ababa"  "Antananarivo"

attr(,"na.action")
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
attr(,"class")
[1] "omit"

afrotropic_cities_nmds = community_nmds(communities %>% filter(city_id %in% afrotropic_cities_community_data$city_id)) 

Run 0 stress 0.00009014786 
Run 1 stress 0.00009691901 
... Procrustes: rmse 0.0001926718  max resid 0.0004158943 
... Similar to previous best
Run 2 stress 0.0003619111 
... Procrustes: rmse 0.001981407  max resid 0.003127616 
... Similar to previous best
Run 3 stress 0.00009161166 
... Procrustes: rmse 0.0002339261  max resid 0.0003170711 
... Similar to previous best
Run 4 stress 0.001912232 
Run 5 stress 0.001215462 
Run 6 stress 0.0006860128 
Run 7 stress 0.00009535671 
... Procrustes: rmse 0.00001612274  max resid 0.00002766371 
... Similar to previous best
Run 8 stress 0.00009163152 
... Procrustes: rmse 0.0002274177  max resid 0.0002943484 
... Similar to previous best
Run 9 stress 0.001472905 
Run 10 stress 0.3017338 
Run 11 stress 0.00009434023 
... Procrustes: rmse 0.0001876868  max resid 0.0004059187 
... Similar to previous best
Run 12 stress 0.3017338 
Run 13 stress 0.00009964034 
... Procrustes: rmse 0.0004865478  max resid 0.001033375 
... Similar to previous best
Run 14 stress 0.002511588 
Run 15 stress 0.00009880586 
... Procrustes: rmse 0.0001946382  max resid 0.0004211662 
... Similar to previous best
Run 16 stress 0.00009893091 
... Procrustes: rmse 0.0001949213  max resid 0.0004193287 
... Similar to previous best
Run 17 stress 0.000468451 
... Procrustes: rmse 0.002599256  max resid 0.004076805 
... Similar to previous best
Run 18 stress 0.0009880508 
Run 19 stress 0.00009746021 
... Procrustes: rmse 0.0001852616  max resid 0.000383588 
... Similar to previous best
Run 20 stress 0.00009997242 
... Procrustes: rmse 0.000481814  max resid 0.001031423 
... Similar to previous best
*** Best solution repeated 12 times

Warning: stress is (nearly) zero: you may have insufficient data

afrotropic_cities_nmds

scree_plot(afrotropic_cities_nmds)

afrotropic_cities = cluster_cities(city_nmds = afrotropic_cities_nmds, cities_community_data = afrotropic_cities_community_data, centers = 4)

Joining with `by = join_by(city_id)`

plot_nmds_clusters(afrotropic_cities)

afrotropic_biomes = resolve[resolve$REALM == 'Afrotropic',c('REALM')]
 
ggplot() + 
  geom_sf(data = afrotropic_biomes, aes(geometry = geometry)) + 
  geom_sf(data = afrotropic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_clusters.jpg'))

Saving 7.29 x 4.51 in image

Afrotropic Cluster 1

afrotropic_cities %>% filter(cluster == 1) %>% plot_city_cluster('Afrotropic cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster1.jpg'))

Saving 7.29 x 4.51 in image

Afrotropic Cluster 2

afrotropic_cities %>% filter(cluster == 2) %>% plot_city_cluster('Afrotropic cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster2.jpg'))

Saving 7.29 x 4.51 in image

Afrotropic Cluster 3

afrotropic_cities %>% filter(cluster == 3) %>% plot_city_cluster('Afrotropic cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster3.jpg'))

Saving 7.29 x 4.51 in image

Afrotropic Cluster 4

afrotropic_cities %>% filter(cluster == 4) %>% plot_city_cluster('Afrotropic cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster4.jpg'))

Saving 7.29 x 4.51 in image

Indomalayan

indomalayan_cities_community_data = community_data_metrics %>% filter(core_realm == 'Indomalayan')
indomalayan_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()

$city_name
  [1] "Srinagar"            "Jamnagar"            "Jammu"               "Rajkot"              "Bikaner"             "Jodhpur"            
  [7] "Jalandhar"           "Ahmedabad"           "Bhavnagar"           "Ludhiana"            "Anand"               "Udaipur"            
 [13] "Surat"               "Vadodara"            "Ajmer"               "Chandigarh"          "Vasai-Virar"         "Mumbai"             
 [19] "Jaipur"              "Delhi [New Delhi]"   "Nashik"              "Dehradun"            "Kota"                "Pune"               
 [25] "Haridwar"            "Dhule"               "Ujjain"              "Indore"              "Ahmadnagar"          "Kolhapur"           
 [31] "Jalgaon"             "Agra"                "Aurangabad"          "Sangli"              "Belagavi"            "Gwalior"            
 [37] "Budaun"              "Bareilly"            "Dharwad"             "Bhopal"              "Bhind"               "Mangaluru"          
 [43] "Solapur"             "Vijayapura"          "Akola"               "Latur"               "Kannur"              "Davanagere"         
 [49] "Thalassery"          "Amravati"            "Kalaburagi"          "Kozhikode"           "Guruvayur"           "Malappuram"         
 [55] "Lucknow"             "Thrissur"            "Mysuru"              "Kochi"               "Alappuzha"           "Nagpur"             
 [61] "Kollam"              "Jabalpur"            "Ettumanoor"          "Hyderabad"           "Coimbatore"          "Bengaluru"          
 [67] "Thiruvananthapuram"  "Tiruppur"            "Faizabad"            "Erode"               "Prayagraj"           "Pratapgarh"         
 [73] "Salem"               "Dindigul"            "Madurai"             "Tiruchirappalli"     "Durg"                "Vellore"            
 [79] "Tirupati"            "Raipur"              "Bilaspur"            "Vijayawada"          "Puducherry"          "Chennai"            
 [85] "Kathmandu"           "Colombo"             "Rajamahendravaram"   "Patna"               "Kandy"               "Bihar Sharif"       
 [91] "Visakhapatnam"       "Ranchi"              "Brahmapur"           "Jamshedpur"          "Darjeeling"          "Siliguri"           
 [97] "Cuttack"             "Bhubaneshwar"        "Jalpaiguri"          "Berhampore"          "Kolkata"             "Krishnanagar"       
[103] "Guwahati [Dispur]"   "Agartala"            "Silchar"             "Dimapur"             "Bangkok"             "George Town"        
[109] "Kuala Lumpur"        "Phnom Penh"          "Singapore"           "Hong Kong"           "Sha Tin"             "Hsinchu"            
[115] "Taichung"            "New Taipei [Taipei]" "Tainan"              "Denpasar"            "Kaohsiung"           "Kota Kinabalu"      

attr(,"na.action")
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
attr(,"class")
[1] "omit"

indomalayan_cities_nmds = community_nmds(communities %>% filter(city_id %in% indomalayan_cities_community_data$city_id)) 

Run 0 stress 0.1190668 
Run 1 stress 0.119963 
Run 2 stress 0.1616726 
Run 3 stress 0.1242016 
Run 4 stress 0.1582327 
Run 5 stress 0.1183841 
... New best solution
... Procrustes: rmse 0.006657586  max resid 0.04580956 
Run 6 stress 0.1219562 
Run 7 stress 0.1153501 
... New best solution
... Procrustes: rmse 0.02592979  max resid 0.2352578 
Run 8 stress 0.1509741 
Run 9 stress 0.155088 
Run 10 stress 0.128717 
Run 11 stress 0.1161599 
Run 12 stress 0.1526359 
Run 13 stress 0.1528153 
Run 14 stress 0.1373095 
Run 15 stress 0.1168773 
Run 16 stress 0.152169 
Run 17 stress 0.1296541 
Run 18 stress 0.1256725 
Run 19 stress 0.1526831 
Run 20 stress 0.1364842 
Run 21 stress 0.1153499 
... New best solution
... Procrustes: rmse 0.00005138461  max resid 0.0003304756 
... Similar to previous best
*** Best solution repeated 1 times

indomalayan_cities_nmds

scree_plot(indomalayan_cities_nmds)

indomalayan_cities = cluster_cities(city_nmds = indomalayan_cities_nmds, cities_community_data = indomalayan_cities_community_data, centers = 6)

Joining with `by = join_by(city_id)`

plot_nmds_clusters(indomalayan_cities)

indomalayan_biomes = resolve[resolve$REALM == 'Indomalayan',c('REALM')]
 
ggplot() + 
  geom_sf(data = indomalayan_biomes, aes(geometry = geometry)) + 
  geom_sf(data = indomalayan_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_clusters.jpg'))

Saving 7.29 x 4.51 in image

Indomalayan Cluster 1

indomalayan_cities %>% filter(cluster == 1) %>% plot_city_cluster('Indomalayan cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster1.jpg'))

Saving 7.29 x 4.51 in image

Indomalayan Cluster 2

indomalayan_cities %>% filter(cluster == 2) %>% plot_city_cluster('Indomalayan cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster2.jpg'))

Saving 7.29 x 4.51 in image

Indomalayan Cluster 3

indomalayan_cities %>% filter(cluster == 3) %>% plot_city_cluster('Indomalayan cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster3.jpg'))

Saving 7.29 x 4.51 in image

Indomalayan Cluster 4

indomalayan_cities %>% filter(cluster == 4) %>% plot_city_cluster('Indomalayan cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster4.jpg'))

Saving 7.29 x 4.51 in image

Indomalayan Cluster 5

indomalayan_cities %>% filter(cluster == 5) %>% plot_city_cluster('Indomalayan cluster 5')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster5.jpg'))

Saving 7.29 x 4.51 in image

Indomalayan Cluster 6

indomalayan_cities %>% filter(cluster == 6) %>% plot_city_cluster('Indomalayan cluster 6')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster6.jpg'))

Saving 7.29 x 4.51 in image